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We propose a method to study lattice QCD at non-vanishing temperature (T) 
and chemical potential (/i). We use n^=2+l dynamical staggered quarks with 
i semi-realistic masses on Lt=4 lattices. The critical endpoint (E) of QCD on the 

Re(/x)-T plane is located. We calculate the pressure (p), the energy density (e) 
and the baryon density (ng) of QCD at non- vanishing T and fi. 
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1. Introduction 

QCD at finite T and/or fi is of fundamental importance, since it describes 
physics relevant in the early universe, in neutron stars and in heavy ion 
collisions. Extensive experimental work has been done with heavy ion col- 
lisions at CERN and Brookhaven to explore the fi-T phase boundary (cf. 
1 ). Note, that past, present and future heavy ion experiments with always 
higher and higher energies produce states closer and closer to the T axis of 
the (jL-T diagram. It is a long-standing question, whether a critical point ex- 
ists on the /j-T plane, and particularly how to tell its location theoretically 

2 

Universal arguments 3 and lattice results 4 indicate that at fj,=0 the real 
world probably has a crossover. Arguments based on a variety of models 
(see e.g. 6 - 7 > 2 ) predict a first order finite T phase transition at large /z. 
Combining the [i — and large /i informations suggests that the phase 
diagram features a critical endpoint E (with chemical potential he and 
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temperature Tg), at which the line of first order phase transitions (/i > He 
and T <Te) ends 2 . At E the phase transition is of second order and long 
wavelength fluctuations appear, which results in (see e.g. 8 ) consequences, 
similar to critical opalescence. Passing close enough to (he,Te) one expects 
simultaneous signatures which exhibit non-monotonic dependence on the 
control parameters 9 , since one can miss the critical point on either side. 

The location of E is an unambiguous, non-pcrturbative prediction of 
QCD. No ab initio, lattice QCD study based was done to locate E. Crude 
models with m s — oo were used (e.g. 2 ) suggesting that he ~ 700 MeV, 
which should be smaller for finite m s . The goal of our work is to propose 
a new method to study lattice QCD at finite /i and apply it to locate the 
cndpoint. We use full QCD with dynamical nf— 2+1 staggered quarks. 

QCD at finite \i can be given on the lattice 10 ; however, standard Monte- 
Carlo techniques fail. At Re(i<)^0 the determinant of the Euclidean Dirac 
operator is complex, which spoils any importance sampling method. 

Several suggestions were studied in detail to solve the problem. 

For small gauge coupling an attractive approach is the "Glasgow 
method" 11 in which the partition function is expanded in powers of 
exp(/u/T) by using an ensemble of configurations weighted by the ti=0 ac- 
tion. After collecting more than 20 million configurations only unphysical 
results were obtained: a premature onset transition. The reason is that the 
H=0 ensemble does not overlap sufficiently with the states of interest. We 
show how to handle this problem for small h values. 

At imaginary fi the measure remains positive and standard Monte Carlo 
techniques apply. One can also use the fact that the partition function away 
from the transition line should be an analytic function of Hi and the fit for 
imaginary \x values could be analytically continued to real values of ji 12 > 13 . 

Due to the renewed interest several promising ideas appeared in the last 
years (without giving a complete list see e.g. 14 . 15 . 16 . 17 . 18 ). 

We propose a method to reduce the overlap problem and determine 
the phase diagram in the /x-T plane (for details see 19 ). The idea is to 
produce an ensemble of QCD configurations at /x=0 and at the transition 
temperature T c . Then we determine the Boltzmann weights 20 of these con- 
figurations at fj, ^ and at T lowered to the transition temperatures at this 
non-vanishing it. Since transition configurations are reweighted to transi- 
tion ones a much better overlap can be observed than by reweighting pure 
hadronic configurations to transition ones 11 . Since the original ensemble 
is collected at /x=0 we do not expect it to be able to describe the physics of 
the large \x region with e.g. exotic colour superconductivity. Fortunately, 
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Figure 1. Schematic difference between the present and the Glasgow methods. 

the typical /x values at present heavy ion accelerators are smaller than the 
region we cover. 

We illustrate the applicability of the method and locate the critical point 
of QCD. (Multi-dimensional reweighting was successful for determining the 
endpoint of the hot electroweak plasma 22 e.g. on 4D lattices.) Further- 
more we give the equation of state of the QCD plasma at non-vanishing 
temperature and chemical potential. 

2. Overlap improving multi-parameter reweighting 

Let us study a generic system of fermions tp and bosons <j>, where the fermion 
Lagrange density is {pM((j))tp. Integrating over the Grassmann fields we get: 



where a is the parameter set of the Lagrangian. In the case of staggered 
QCD a consists of (3, m q and it. For some choice of the parameters a=a n 
importance sampling can be done (e.g. for Re(/i)=0). Rewriting eq. (1) 



The curly bracket is measured on each independent configuration and is 
interpreted as a weight factor {w(f3, //, m, £/)}. The rest is treated as the 
integration measure. Changing only one parameter of the ensemble gener- 
ated at ao provides an accurate value for some observables only for high 




(1) 
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statistics. This is ensured by rare fluctuations as the mismatched measure 
occasionally sampled the regions where the integrand is large. This is the 
overlap problem. Having several parameters the set «o can be adjusted to 
get a better overlap than obtained by varying only one parameter. 

The basic idea of the method as applied to dynamical QCD can be 
summarized as follows. We study the system at Re(/z)=0 around its transi- 
tion point. Using a Glasgow- type technique we calculate the determinants 
for each configuration for a set of /i, which, similarly to the Ferrenberg- 
Swendsen method 20 , can be used for reweighting. The average plaquette 
values can be used to perform an additional reweighting in (3. Since transi- 
tion configurations were rewcightcd to transition ones a much better overlap 
can be observed than by reweighting pure hadronic configurations to transi- 
tion ones as done by the Glasgow- type techniques. The differences between 
the two methods are shown in Figure 1. (Note, that reweighting techniques 
have very broad applicability. E.g. recently it was possible to determine 
the topological susceptibility with overlap fermions 21 .) 

3. The endpoint of n f — 2 + 1 QCD 

In QCD with rif staggered quarks one changes the determinants to their 
rtf/A power in our two equations. Importance sampling works at some (3 
and at Ke(/j,)—0. Since detM is complex an additional problem arises, one 
should choose among the possible Ricmann-sheets of the fractional power 
in eq. (2). This can be done by using 19 the fact that at fi=^ w the ratio of 
the determinants is 1 and it should be a continuous function of [i. 

In the following we keep n real and look for the zeros of Z for complex /?. 
At a first order phase transition the free energy oc log Z{(3) is non-analytic. 
A phase transition appears only in the oo limit, but not in a finite 
V. Nevertheless, Z has zeros at finite V, which are at complex parameters 
(e.g. (3). For a system with first order transition these zeros approach the 
real axis as oo by a 1/V scaling. This oo limit generates the non- 
analyticity of the free energy. For a system with crossover Z is analytic, 
and the zeros do not approach the real axis as oo. 

At T^O we used L t —A, L s =4,6,8 lattices. T=0 runs were done on 10 3 - 
16 lattices. m u ,d= 0.025 and m s =0.2 were our bare quark masses. At T ^ 
we determined the complex valued Lee- Yang zeros 23 , (3$, for different V-s 
as a function of /j,. Their oo limit was given by a (3o(V) — f3^° + (/V 
extrapolation. We used 14000, 3600 and 840 configurations on L s =4,6 and 
8 lattices, respectively. For small \x values Im(/3g°) is inconsistent with zero, 
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Figure 2. The T-/i diagram. Direct results are given with errorbars. The dotted line 
shows the crossover, the solid line the first order transition. The box gives the uncer- 
tainties of the endpoint. 

and predicts a crossover. Increasing /z, the value of Im(/3g°) decreases. Thus 
the transition becomes consistent with a first order phase transition. Our 
primary result is fiend = 0.375(20) in lattice units. 

To set the physical scale we used an average of i? , m p and y/a. Includ- 
ing systematics due to finite V we have (i?o • m Tr) = 0.73(6), which is at 
least twice, m u .d is at least four times as large as the physical values. 

Figure 2 shows the phase diagram in physical units, thus T as a function 
of hb, the baryonic chemical potential (which is three times larger than 
the quark chemical potential). The endpoint is at Te = 160 ± 3.5 MeV, 
He = 725 ± 35 MeV. At fi B =0 wc obtained T c = 172 ± 3 MeV. 

Using a Taylor expansion around /x=0, T^O for small \i can be used 
to determine the curvature of the phase diagram and to calculate thermal 
properties 24 . A different method, analytic continuation from imaginary /i, 
confirmed also the results of 19 on the pi-T diagram 12 ' 13 . 

4. Equation of state at non-vanishing T and fi 

The equation of state (EOS) at fi^O is essential to describe the quark 
gluon plasma (QGP) formation at heavy ion collider experiments. Results 
are only available for fi=0 (e.g. 25 > 26 > 27 ) at T^0. 

We use 4 -N% lattices at T^0 with =8, 10, 12 for rcweighting and we 
extrapolate to V^oo using the available volumes (V). At T—0 lattices of 
24 • 14 3 are taken for vacuum subtraction and to connect lattice parameters 
to physical quantities. 14 different (3 values are used, which correspond to 
T/T c — 0.8, . . . , 3. Our T=0 simulations provided Rq and a. The lattice 



7, 2008 9:51 WSPC/Trim Size: 9in x 6in for Proceedings 



heidelberg'02b 



6 




Figure 3. The best weight lines on the fi-[3 plane. In the middle we indicate the 
transition line. Its first dotted part is the crossover region. The blob represents the 
critical endpoint, after which the transition is of first order. The integration paths used 
to calculate p are shown by the arrows along the (3 axis and the best weight lines. 




Figure 4. p normalised by T 4 as a function of T/T c at fi = (to help the continuum 
interpretation the raw lattice result is multiplied by £^=0.446). 



spacing at T c is «0. 25-0. 30 fm. We use 2+1 flavours of dynamical stag- 
gered quarks. While varying /3 (thus T) we keep the physical quark masses 
approx. constant (the pion to rho mass ratio is m w /m p «0.66). 

The determination of the equation of state at /i 7^0 needs several ob- 
servables, O, at fi^O. This is obtained by using the weights of eq. (2) 

0(/3,ii,m) = w 7-5 r+7i ■ 3 

p can be obtained from the partition function as p=T ■ d\ogZ/dV 
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Figure 5. Ap = p(jj, ^ 0, T) - p(n = 0, T) normalised by T 4 as a function of T/T c for 
AtB =100, 210, 330, 410 MeV and 530 MeV (from bottom to top). To help the continuum 
interpretation the raw lattice result is multiplied by c M =0.446. 



which can be written as p=(T/V)TogZ for large homogeneous systems. 
On the lattice we can only determine the derivatives of log Z with respect 
to the parameters of the action (0, m, fx). Using the notation (0((3, fi, m)) = 
0(/3, \i, m) T _£ — 0(/3, \i = 0, m) T=Q . p can be written as an integral 28 : 



P_ 



1 



T 3 V . 
dQogZ) 



d((3,m,fj.) 
d(\ogZ) 



(4) 



dp I ' \ dm I ' \ d(i 
The integral is by definition independent of the integration path. The 
chosen integration paths are shown in Fig 3. 

The energy density can be written as e = {T 2 /V)-d{\og Z) / dT+(pT /V)- 
d(\og Z)jd[i. By changing the lattice spacing T and V are simultaneously 
varied. The special combination e— 3p contains only derivatives with respect 
to a and fj,: 



3 P 



d\og(Z) 



T 4 



T 3 V da 



V d(\ogZ) 



T 3 V dfi 



(5) 



The quark number density is n — (T/V) ■ d log(Z) / '<9/i which can be mea- 
sured directly or obtained from p (baryon density is ns=n/3 and baryonic 
chemical potential is /zb=3/u). 

We present direct lattice results on = 0,T), Ap(fi,T) = p(p, ^ 
0,T)-p(n = 0,T), e(/x, T)-3p(/i, T) and n B {fi,T). Additional overall fac- 
tors were used to help the phenomenological interpretation. 
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Figure 6. (e - 3p)/T 4 at /i fl =0, 210, 410 MeV and 530 MeV versus T/T c (from bottom 
to top). To help the continuum interpretation the raw lattice result is multiplied by 
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Figure 7. n s /T 3 versus T/T c for the same /ig values as in Fig. 3 (from bottom to top), 
(to help the continuum interpretation the raw lattice result is multiplied by c,j=0.446). 
As a reference value the line starting in the left upper corner indicates the nuclear density. 



Fig. 4 shows p at /j,=0. In Fig. 5 we show Ap/T 4 for different \i values. 
Fig. 6 shows e-3p normalised by T 4 , which tends to zero for large T. Fig. 
7 gives the baryonic density as a function of T/T c for different fjrS. The 
densities can exceed the nuclear density by up to an order of magnitude. 

An important finding concerns the applicability of our reweighting 
method: the maximal fj, scales with the volume as /i max • a ~ (N t ■ Ng)~° . 
If this behaviour persists, one could -in principle- approach the true con- 
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tinuum limit (a ~ 1/N t ~ (iV 4 • N^)^ - 25 , thus ^ max «const.). 
5. Conclusion 

Wc proposed a method an overlap improving multi-parameter reweighting 
technique- to numerically study non-zero \x and determine the phase dia- 
gram in the T-\i plane. Our method is applicable to any number of Wilson 
or staggered quarks. 

We studied the fi-T phase diagram of QCD with dynamical n/=2+l 
quarks. Using our method we obtained Te«160 MeV and /xe«700 MeV 
for the endpoint. Though \ie is too large to be studied at RHIC or LHC, 
the endpoint would probably move closer to the /i=0 axis when the quark 
masses get reduced. 

The equation of state was determined on the temperature versus chem- 
ical potential plane. According to our results the applicability range of 
the overlap improving multi-parameter reweighting method for the quark 
chemical potential can be summarized as n<T. 

Clearly, more work is needed to get the final values by extrapolating in 
the R-algorithm and to the thermodynamic, chiral and continuum limits. 
The details of the presented results can be found in 19 ' 29 . 
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